Monte Carlo computation of pair correlations in excited nuclei 
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We present a novel quantum Monte Carlo method based on a path integral in Fock space, which 
allows to compute finite-temperature properties of a many-body nuclear system with a monopole 
pairing interaction in the canonical ensemble. It enables an exact calculation of the thermodynamic 
variables such as the internal energy, the entropy, or the specific heat, from the measured moments 
of the number of hops in a path of nuclear configurations. Monte Carlo calculations for a single- 
shell (/iii/2)^ model are consistent with an exact calculation from the many-body spectrum in the 
seniority model. 
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In the last dorado, many advances in the solution of the 
nuclear many-body problem have become possible due to 
the development of quantum Monte Carlo (MC) meth- 
ods. Standard Green's Fimction Monte Carlo (GFMC) 
methods have been used to perform a microscopic calcu- 
lation of the ground state of nuclei up to A = 6 with real- 
istic interactions [1.2]. Recently, a path-integral cjuantum 
Monte Carlo method has also been devised to explore 
the microscopic properties of nuclei at finite tempera- 
ture [3]. This method, called Shell Model Monte Carlo 
(SMMC), uses the Hubbard-Stratonovich transformation 
of the many-body propagator to express observables as 
path integrals of one-body propagators in fluctuating 
auxiliary fields. It has been used successfully to treat 
shell-model spaces larger than those accessible to conven- 
tional methods [4-6]. In this Letter, we explore an alter- 
native quantum MC method, based on a path integral in 
Fock space, that allows the calculation of thermodynamic 
properties of the many-body system. The observables are 
obtained through the MC computation of path integrals 
in Fock space, rather than in a space of auxiliary field 
variaV)les as in SMMC, or in a spatial (and spin/isospin) 
coordinate space as in GFMC. One main advantage is 
that it yields a proper description of excited nuclei in 
terms of a canonical ensemble, without having to be con- 
cerned with particle-number projection. The method 
permits considering large model spaces, since it avoids 
an explicit enumeration of the many-body states. The 
results are in principle exact, i.e. they are subject only 
to controllable sampling errors. The time-discretization 
error, inherent to traditional path-integral MC methods 
(see e.g. [4,5]), is absent here due to the discrete nature 
of the states in Fock space. 

An important motivation for an exact treatment of 
thermal properties of nuclei comes from the description 
of low- energy reactions, which depend crucially on level 
densities. In situations where the levels cannot be mea- 
sured, as for example in the very neutron-rich nuclei 
involved in r-process nucleosynthesis, one still relies on 
semi-plienomenological level densities to calculate cross- 
sections [7]. In this Letter, we study a specific nuclear 
interaction, the pairing force, which exhibits the simplic- 
ity and efficiency of the proposed MC technique. It is 
well known that the finite-temperature BCS theory pre- 
dicts a spurious phase transition due to particle-number 
fiuctuations: this transition occurs at a temperature cor- 
responding to an excitation energy of interest in many 
cases. Therefore, even if more refined methods have been 
devised such as the static-path approximation [8], an ex- 
act treatment is of a great help. The present work offers 
the perspective that, using an adecjuate procedure to cir- 
cumvent the "sign problem" generic to fermionic systems 
[9], more general interactions could also be described by 
our path-integral MC method in Fock space; this is under 
investigation. 

Consider a 7V-nucleon system with a monopole pairing 



interaction, described by the Hamiltonian 

Q n 

k = l k,k' = l 
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where k and —k are time-reversed conjugate orbitals of 
energy f j.. Our MC simulation is based on a random walk 
in Fock space, where the two-body pairing interaction can 
be seen as a pair hopping term [10]: a pair of nucleons in 
orbits (fc, —k) hops to (fc'. —k') with an amplitude Gfcj.'. 
The set of basis states of our Fock space, i.e. configura- 
tions of A'^ nucleons into 2fl orbitals, are noted C in the 
following. Because of the particular form of H. the hop- 
ping term does not act on unpaired nucleons; therefore, 
in this C-representation, one is considering a "mixed" 
system, where unpaired nucleons behave like "classical" 
particles, while pairs are "boson-like" particles (with a 
hard-core interaction preventing two pairs from occupy- 
ing the same orbit pair). The thermodynamic properties 
of this system are determined from a path-integral MC 
computation of the density matrix. The decomposition 
of that matrix in C-space. p{C.C') = (C|e~''^|C"), sat- 
isfies the so-called "sign rule", i.e. all matrix elements 
are non-negative, because the pairing interaction is at- 
tractive (Gk^k' > 0). As a consequence, this system can 
be efficiently simulated by quantum MC in C-space. In 
contrast, when treating more general interactions, paths 
of positive and negative sign that tend to cancel each 
other would contribute to the path integral in Fock space, 
thereby requiring to devise a specific solution to that 
"sign proV)lem" for efficiently simulating the system. 

In order to treat the operator e~^^ in C- 
representation, we work in interaction picture, writing 
H = Ho + H\ where Hg and Hj are the diagonal and 
off-diagonal part of H, respectively. We start with the 
Green's function G(C, C";/) in imaginary time (h = 1), 
solution of the integral equation 

G{t) = Go{t)- f dt'Go{t')HiG{t-t') (2) 
Jo 

where Go = e^*^°, diagonal in C-space. plays the role 
of a "free-particle" propagator. The pair hopping term 
-f^i = Y^kj^k' Gk,k'a\,a^-k''^-k'^k can be seen as an "in- 
teraction" Hamiltonian. Applying Erj. (2) recursively, as 
usual, yields an expansion of the full propagator in terms 
of number of pair hops 

G{t) = Go{t)- f dtiGo{ti)HiGo{t - ti)+ (3) 
Jo 

[ AtiAt2Go(ti)HiGo{t2 - ti)HiGo(t -t^) 

This expansion suggests one can perform integration over 
the paths in C-space with a MC simulation simply by 
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sampling the times f^. f-j. - ■ ■ at which pair hops occur, 
and the configurations C in each inter-hop interval. 

Consider now the MC calculation of thermodynamic 
averages in the canonical ensemble. The thermal expec- 
tation value of an observable O at inverse temperature p, 
{0)p = Tr(Oe^^^)/Tr(e^^^), is derived by a stochastic 
calculation of the partition function Z(f3) = Tr(e~^^ ). 
Using the expansion of the many-body Green's function 
G(C\ C"; t) in C-space, Z(j3) can be written as a path 
integral over chains of configurations, noted C(f). peri- 
odic in time (of period /3), having pair hops at times 
ti, ^2, ■ ■ ■ ^-D- Here D = D[C'(t)] is the number of hops 
[11] in the chain C'{t), and it can be any non-negative 
integer excluding one. It is the periodic boundary con- 
dition which forbids D = I, and we will see below that 
it has important consequences. Note that, due to par- 
ticle indistinguishability, periodic chains which make an 
overall pair permutation in one period have also to be 
considered. Thus, in order to calculate thermodynamic 
averages by MC, one has to sample the whole ensemble 
of chains of configurations {C'{t)} with variable D ^ 1 
by use of the Metropolis algorithm [12]. That algorithm 
permits to sample a chain C'{t) with a probability pro- 
portional to its weight iy[C(t)] in Z{J3), given from Eq. 
(3) by 

D 

W[C{t)]dhdt2 ■■■dtD = exp {A[C(t)]} n Gk,k'{i)dU 

(4) 

Here A[C(i)] = - E^'' (C (t))dt is the action in imagi- 
nary time corresponding to the one-body energy E^'' of 
the configuration chain, and Gk,k'{i) is the amplitude as- 
sociated with the i-th hop along the chain. Since the 
propagator does not hop unpaired nucleons, the latter 
stay on the same orbit along the chain; their contribu- 
tion to E^^ is constant in time, simply giving an additive 
constant to the action. 

Now, making use of the symmetry in imaginary time 
- all possible times t at which the observable O is in- 
serted in the chain C(i) when calculating {O) g are equiv- 
alent - one can write thermal expectation values as av- 
erages over paths and over insertion times r. This gives 
H\C{t)] = E[C'{t)] — D/ff as an estimator for the energy, 
with E[C{t)] = E^''{C{t))dt being the one-body 

energy averaged along the chain C'(t). An estimate of the 
internal energy is then obtained by 

U = {H)p - So ^ mem) at) - E, (5) 

where Eq stands for the ground state energy (extracted 
from the /3 — ^ oo limit) and {{■)) denotes a MC aver- 
age over the sampled chains C{t). Thus, U is a func- 
tion of the average number of hops along the chain 
{{D)) (any additional hop gives on average the con- 
tribution to U). The same reasoning yields an 



estimate of the pair occupation probability in the k- 
th orbital, (n,,),3 ~ {{nk\C{t)\))c(t), with nk\C{t)\ = 
nk{C{f))dt. Thus, the occupation of the fc-th or- 
bital, averaged along each chain C(f) and averaged over 
the sampled chains, simply yields the exact quantum ex- 
pectation value of Uk at inverse temperature j3. One can 
also obtain an estimator of the specific heat from the 
variance of the energy, C = /3^Var{i?}, which is given by 

Cc^p'Yar{H[C{t)]}at)-{{D))c(t) (6) 

Thus, the specific heat depends on the variance of the 
number of hops along the chain, and on the correla- 
tion between the number of hops D and the averaged 
one-body energy of a chain E[C(t)]. Finally, the en- 
tropy S in the canonical enseml)le can be estimated 
by use of MC (see below). Then, the level density 
can be calculated via the saddle-point approximation as 
p{U) = :3c^ / {'InCy using the MC estimates for S and 
C. This provides a natural extension to the MC approach 
to level densities for non-interacting nucleons suggested 
in [13]. 

A general difficulty related to the Metropolis method 
is to find an efficient algorithm for generating trial moves 
(i.e., for creating a new C'(t) from the current one), 
which has a non-negligible acceptance rate and a rea- 
sonably short auto-correlation time. The detail of a 
general algorithm for sampling the paths C(t) will be 
reported elsewhere. For illustrative purposes, we fo- 
cus in this Letter to the case of a single-shell model 
(ej, = 0) with a constant pairing strength G. This 
model, called the seniority model [14], can been solved 
exactly using quasi-spin formalism, allowing us to demon- 
strate the interest of the MC procedure. Defining the 
seniority quantum number s (s = 0,2, - ■■ N for even 
only considered here), the A^-body energy levels are 
E{N,s) = -G{N - s)(2Sl - s - N + 2)/4 with a de- 
generacy d(s) = (^"2) ~ is/2--i)- ^'^^^ '''(^) = 1- This 
model is simpler to compute by MC than a general case 
because £'^'((7) = —NGj'l is independent of C. so that 
the weight of any chain having D hops scales as G". 
Also, integration over (time-ordered) intermediate times 
t\. f'l ' ' ' to f'fin be made explicitly (without rerjuiring MC 
calculation), giving the factor / D\ in the weigths, so 
that it is enough to sample periodic sequences of C's of 
variable length D (D = 0. 2.3. •••)• In those sequences 
{Ci, C2, ■ ■ ■ C_d}, Ci and Ci+i must differ by exactly one 
hop. 

As an illustration, we consider a half- filled /111/2 shell 
with a constant pairing strength G = 1 MeV. We plot in 
Fig. 1 the thermal expectation value of energy as a func- 
tion of inverse temperature. The energy in the seniority 
model is estimated by MC from -NG /2- p''^ {{D)). The 
sfjuares are the result of MC simulations using a sample 
of 10^ sequences, while the solid line was calculated nu- 
merically from the exact spectrum in the seniority model. 
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FIG. 1. Temperature dependence of the expectation value 
of the many-body energy, {H)^, in a single-shell (/111/2)* 
model with a constant pairing strength G = 1 MeV. The 
squares represent Monte Carlo results at each inverse tem- 
perature /3 = T^^. The solid curve corresponds to an 
exact numerical calculation from the spectrum in the se- 
niority model. In the inset, we show the average number 
of hops {{D)) derived from the MC simulations (squares) 
as a function of the chain length /3. The dashed curve, 
corresponding to the asj-mptotic (/3 — > 00) behaviour 
{{D)) = fi = l3GN (2n - N)/A, illustrates the hindrance due 
to periodic boundary conditions. 
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The agreement is very good. The inset in Fig. 1 shows 
the corresponding average number of hops {{D)) in the 
chains sampled by MC (squares) as a function of the 
chain length '3. For long chains. {{D)) increases linearly 
with the chain length, with the asymptotic ( 3 — > 00) be- 
haviour {{D)) = l3GN{2ft - N)/4 (dashed line). This 
is expected since it was shown in rof. [10] that, at zero 
temperature, the number of hops [15] per imaginary-time 
unit At is drawn from a Poisson distribution of param- 
eter u = AtGN(2Q — N)jA. For short chains, however, 
{{D)) increases quadratically, as j3'^G'^N{2Vt — N)j'i, and 
there is a depletion in the average number of hops. That 
depletion, which reflects the hindrance due to periodic 
boundary conditions (D = 1 is excluded), implies the 
non-zero specific heat in this model as shown below. 

In Fig. 2, we plot the entropy 5* as a function of tem- 
perature T = f3~^. The Monte Carlo calculation of S 
in the canonical ensemble (squares) uses relation S = 
\nZ{p) + 3U , where In Z(3) is obtained by numerically 
integrating U from to /3 (remember din Zj dp = —U), 
using S{0) = In (j^j) ■ T^^^ agreement with the exact 
numerical calculation of S (solid line) derived from the 



FIG. 2. Temperature dependence of the entropy, S, calcu- 
lated from S = lnZ{0)-\- /SU in the same model. The Monte 
Carlo calculations correspond to squares, while the exact cal- 
culation is shown as a solid curve. In the inset, we show the 
corresponding results for the specific heat C, calculated by a 
finite difference approximation to dU/dT. The observed peak 
is a clear signature of the vanishing of pair correlations at a 
temperature around T ~ 2 MeV. 
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spectrum is evident. The inset in Fig. 2 shows the tem- 
perature dependence of the specific heat C in this model. 
Since E[C(f)] = —NG/2 in the seniority model, the MC 
estimator is simply C ~ Var{73} — {{D)). According to 
ref. [10], for /3 — ^ 00, i.e. when the effect of periodic 
boundary conditions vanishes, the number of hops D in 
the chain of length l3 = MAt is distributed as a Poisson 
variable of parameter jj = Mv = pGN(2fl — N)/4:] then, 
it is obvious that C tends to zero for 3 — > -dc. as expected. 
An independent method to access to the specific heat by 
MC is to make use of C = dU/dT = -f3'^dU/df3. Graphi- 
cally, it implies that the tangent of the curve {{D)) versus 
/3 (see inset of Fig. 1) intercepts the y-axis at — C; there- 
fore, the depletion of that curve compared to the asymp- 
totic curve indicates that G tends to a maximum value 
and then falls to zero. The squares, obtained by MC (by 
numerically estimating the derivative of U), are in good 
agreement with the solid line, calculated from the exact 
spectrum. The MC determination of G from the energy 
variance (not shown in Fig. 2) is also compatible with the 
exact result, but is affected by a much larger statistical 
noise for large /3. The peak in the specific heat around 
T ~ 2 MeV refiects the pairing phase transition; above 
that temperature, the pair correlations vanish. The MC 
simulations also exhibit a collapse in the measured pair- 
ing field (A'l'A) near 2 MeV, with A = a-kctk, which 
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is a clear signature of the pairing transition. 

In conclusion, we have shown that a path-integral 
Monte Carlo method in Fock space is well-suited to 
describing pair correlations in nuclei at finite temper- 
ature. It allows a stochastic calculation of thermody- 
namic properties (internal energy, entropy, specific heat, 
level density) in the canonical ensemble: there is no need 
to project into the number of particles. The scatter- 
ing of pairs in the orbits is simulated straigthforwardly; 
for instance, the pair occupation probabilities are mea- 
sured during the simulation, not calculated. Estimators 
for thermal expectation values of observables have sim- 
ple expressions depending on the number of pair hops D 
in the path; for instance, the internal energy (or specific 
heat) is related to the average (or variance) of D. Results 
for a single-shell (/iii/2)^ model with a constant pairing 
strength of 1 MeV are compared with a numerical calcu- 
lation from the exact spectrum in order to ascertain the 
method. The observed peak in the specific heat is asso- 
ciated with the vanishing of nucleon pair correlations at 
a temperature of ~ 2 MeV. In contrast to conventional 
path-integral MC calculations which need to discretize 
the imaginary time into time slices At (see e.g. [4,5]), 
our MC method uses a continuous-time formalism; this is 
possible because the path integral lies in a discrete config- 
uration space. Therefore, there is no time-discretization 
error (or no need to extrapolate to At — ^ 0) and the 
only residual error in our MC method is the statistical 
noise, which can be made arbitrarily small. The exten- 
sion of these MC calculations to a more general pairing 
interaction, or even to other interactions, should allow an 
exact description of nuclear properties at high excitation 
energies. 
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